Implementation of screened hybrid functionals based on the Yukawa potential within 
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The implementation of screened hybrid functionals into the wien2k code, which is based on the 
LAPW basis set, is reported. The Hartree-Fock exchange energy and potential are screened by means 
of the Yukawa potential as proposed by Bylander and Kleinman [Phys. Rev. B 41, 7868 (1990)] 
for the calculation of the electronic structure of solids with the screened-exchange local density 
approximation. Details of the formalism, which is based on the method of Massidda, Posternak, 
and Baldereschi [Phys. Rev. B 48, 5058 (1993)] for the unscreened Hartree-Fock exchange are given. 
The results for the transition-energy and structural properties of several test cases are presented. 
Results of calculations of the Cu electric-field gradient in CU2O are also presented, and it is shown 
that the hybrid functionals are much more accurate than the standard local-density or generalized 
gradient approximations. 

PACS numbers: 71.15.Ap, 71.15.Dx, 71. 15. Mb 
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INTRODUCTION 



Until now, most Kohn-Sham (KS) density functional 
theory (DFT) 1 -^ calculations on solids have been done 
using either the local density approximation (LDA) or 
the generalized gradient approximation (GGA) for the 
exchange-correlation energy. Calculations were done ex- 
clusively with LDA 2 - until the early 90s, when the GGA 
functional PW91— was proposed and then implemented in 
computer codes for solid-state calculations. A few years 
later, a GGA functional with a simpler analytical form 
than PW91, namely PBE, 4 but giving nearly identical 
results has been proposed and is nowadays the standard 
functional. The successes of the semilocal LDA and GGA 
approximations rely on the fact that the accuracy is usu- 
ally good enough to be useful, in particular for the calcu- 
lation of the geometrical parameters and other quantities 
like the bulk modulus or the phonon spectrum. (See, e.g., 
Refs. [H and[^ for a compilation of lattice constants and 
bulk moduli calculated with various GGA functionals.) 
However, it is known that there are classes of systems, 
e.g., strongly correlated or van der Waals systems, whose 
properties are not described properly by semilocal func- 
tionals already at the qualitative level. 

It is also well known that the KS band gap, defined as 
the conduction band minimum (CBM) minus the valence 
band maximum ( VBM) , obtained from a semilocal func- 
tional is much smaller than the experimental band gap 
(defined as the ionization potential I minus the electron 
affinity A). However, it is important to note that this 
problem, known as the band gap problem, is more gen- 
eral and has its roots in the KS-DFT method itself, and 
actually the KS band gap calculated with the exact mul- 
tiplicative KS potential would differ from I — A by the 
derivative discontinuity A xc of the exchange-correlation 
potential (see Ref. for a review) . Since A xc can be of 
the same order as the KS band gap, the exact KS band 
gap can differ substantially from I — A. 



There are several methods to obtain orbital energies 
which lead to values for CBM — VBM comparable to 
I — A. If one wants to stay inside the true KS frame- 
work (i.e., KS equations with a multiplicative potential), 
exact exchange (EXX) calculations (see, e.g., Refs. |8| 
andlpj) or advanced semilocal potentials 1 ^ can do a good 
job. Alternatively one can use a non-multiplicative po- 
tential, which means to use a method which lies outside 
the KS framework, but belongs to the so-called gener- 
alized KS framework. 11 Most of these methods mix the 
DFT and Hartree-Fock (HF) theories and the best known 
are the LDA+C/^ screened-exchange LDA (sX-LDA)^ 
and hybri d 14 ' 15 methods. The GW method can yield 
very accurate band structures, in particular if it is ap- 
plied self-consistently, but it is a very expensive method 
(see Ref. [H for a review) . 

The LDA+C/ method (see Ref. [l?] for a recent review) 
consists of applying an approximate (but very cheap) 
form of HF only to the electrons which are not well de- 
scribed by semilocal functionals. Typical examples are 
the 3c? or 4/ electrons in strongly correlated systems 
(e.g., transition- metal and rare-earth oxides) that are 
very localized and hence lead to large self-interaction er- 
ror when a semilocal functional is used (this results in too 
small band gaps and magnetic moments). In the sX-LDA 
method, the short-range (SR) part of LDA exchange is re- 
placed by the SR part of HF exchange, where the SR part 
is defined by replacing the bare Coulomb potential by 
the screened Yukawa potential^ into the corresponding- 
expressions for the energy and potential. The sX-LDA 
method has been implemented within the pseudopoten- 
tial plane- wav o 11 ' 13 ' 19 ' 20 and linearized-augmented plane- 
wave 2 ^— basis sets, and it has been shown that sX-LDA 
improves substantially over LDA for the band gap of 
semiconductors and insulators. 

Despite the fact that reports about the implementation 
of the HF method in solid-state codes started to appear 
already in the 70s (see, e.g., Refs. l24l - [27t ). it is only in 
the early 2000s that the first calculations on solids with 
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hybrid methods were reported^ - — which is much later 
than for moleculesj 14 ' 15 In hybrid methods, a certain per- 
centage (between 10% and 50%) of semilocal exchange is 
replaced by HF exchange, while the correlation remains 
purely semilocal. As for molecules, the hybrid function- 
als have shown to lead to (much) better results than 
semilocal functionals for various types of materials and 
properties. In particular, they lead to band structures 
which are usually in good agreement with experiment as 
shown for classical semiconductors and insulators (see, 
e.g., Refs. [28ll29ll3l|) a nd strongly correlated materials 
(see, e.g., Refs. I28ll30l - [34l ) . The most common hybrid 
functionals are B3LYPi^^ and PBEO^l which contain 
20% and 25% of HF exchange, respectively. However, for 
solids the long-range (LR) nature of HF exchange leads to 
technical difficulties. For calculations done in real space, 
the results converge very slowly with respect to the num- 
ber of neighboring unit cells that are taken into account 
for the calculation of HF exchange, while for calculations 
done in reciprocal space, the slow convergence is with re- 
spect to the number of k-points for the integrations in 
the Brillouin zone. 

To reduce this problem of slow convergence, Heyd et 
al. (HSE)— ~— proposed to consider only the SR part 
of HF exchange (as done in sX-LDA), and therefore to 
keep 100% of semilocal LR exchange. This was done by 
splitting the Coulomb operator into SR and LR compo- 
nents by using the error functionj 41 ' 42 Since then, it has 
been shown that the HSE functional, which is based on 
PBE0, leads to very good results for semiconductors and 
insulators 4 ^— including strongly correlated systems 
and several recent papers reporting the implementation 
of HSE have appeared] 44 ' 49 We also mention the onsite 
version of HF exchange proposed by Novak et a/., 50 which 
leads to very cheap calculations, but can be applied only 
to localized electrons. This method has been used in the 
context of hybrid calculations . 51 ' 52 

In the present work, we report the implementation 
of screened hybrid functionals into the WIEN2K code^ 
which is based on the full-potential linearized augmented 
plane-wave plus local orbitals method (abbreviated as 
LAPW in the following) 5 ^— to solve the KS equations. 



As done for the sX-LDA functional^ the HF exchange 
is screened by means of the Yukawa potential in order 
to eliminate the LR HF exchange. The calculation of 
the screened HF exchange is based on the pseudocharge 
method 5 ^ as proposed by Massidda et al. for the un- 
screened HF exchanged In the papers of Asahi et al . 21 ^ 22 
it is mentioned that this mehod was used for the imple- 
mentation of the sX-LDA functional, but only very few 
details are given. At this point we also mention Refs. |59| - 
l62l . in which alternative ways of implementing the HF or 
EXX methods within the LAPW basis set are presented. 

The paper is organized as follows: in Sec. [Til the de- 
tails of the formalism of the unscreened and screened HF 
exchange for the LAPW basis set are given and in Sec. 
IIII1 the implemented screened hybrid functionals are pre- 
sented. In Sec. IIV1 the results for a few test cases and 
CU2O are presented, and in Sec. |V]the summary of the 
work is given. 



II. SCREENED HARTREE-FOCK EXCHANGE 

In this section, the formulas of the screened HF energy 
for the LAPW basis set are given. The formulas are also 
valid (and implemented) for the APW plus local orbitals 
basis seti 56 ' 57 For completeness and to allow comparison, 
the formulas for the unscreened case are also given. For 
the Hamiltonian, only the basic formulas are given. The 
LAPW method will not be described here, but details 
can be found in Refs. I55H571 



A. Energy 

The HF exchange energy per unit cell (of volume f2) 
is given by (all following equations are in Hartree atomic 
units) 
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are the valence- valence (vv), valence-core (vc), and core- 
core (cc) terms, respectively. In Eqs. ([2]) and (J3j> , w° k 
is the product of the k-point weight and the occupation 
number and "0^ k is a spin-cr valence orbital of band index 
n and wave vector k, whose LAPW basis set expansion 
in the interstitial (I) and atomic spheres (S Q ) is given by 
(r a = r — r a , where r Q is the position of nucleus a) 



K 



-n,k+Kn+K 



0k +K (r), 



(5) 



r € I 



1 i(k+K)-r 

0k+ K (r) - S EE^/k+K^«(^)^m(r a ), reS Q 



whcr 



Si.k+K 



(6) 

are the variational coefficients. In the in- 
terstitial, the basis functions </>£ +K are represented by 
plane- waves, while inside the atomic spheres, k+ K are 
linear combinations of products of radial functions u"J 
and spherical harmonics Yg m . The coefficients are 
determined such that the <1>u+k s are continuous across 
the sphere boundaries. For / = 1, 2, and 3, ujl repre- 
sents a radial function evaluated at a linearization energy, 
its energy derivative evaluated at this same energy, and 
a radial function evaluated at another linearization en- 
ergy (e.g., semicore states), respectively. In Eqs. (j3|) and 
©' ^"n ' l m c * s a core or bital which is confined inside 
the atomic sphere S a and where n c , £ c , and m c are the 
principal, azimuthal, and magnetic quantum numbers, 
respectively: 

C c W r ) = <:eSr a )Y £cmc (r a ). (7) 
In Eqs. ([2])-(|4j), v is either the unscreened potential 
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or the Yukawa screened potentia. 



18 



-A r-r 



oo e 



= 4ttA^ i/ (Ar<)fc/(Ar>) 



xy; m (r)r, m (f) 



(9) 



where A is the screening parameter and ii and fc^ are 
spherical modified Bessel functions. 63 Note that the 
spherical harmonics expansion of the screened potential 63 
is simpler than in the case of the error function^ 



1. Valence-valence term 



Following the idea of Massidda et al— the valence- 
valence term [Eq. j2j] is cast into the following form 
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In the interstitial and spheres, p^ kn / k / an( l w nkn'k' are ex ~ 
panded in Fourier and spherical harmonics series, respec- 
tively (from now on, the index a of the position r a from 
the nucleus a is suppressed and we define q = k' k+G): 
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In Eq. (TT3"j) , p^ k „, k , are the Fourier coefficients of the 
periodic part of p£ kn , u , and p^\, is given by 
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with C e ^™^ £m being Gaunt coefficients, 
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and ^r k/ =E K <,k + K^k+k- 

^nkn'k' i s calculated by using Weinert's method for 
solving the Poisson equation.™ (I n Appendix I A 1| a brief 
summary of Weinert's method for the unscreened and 
screened potentials is given.) For the unscreened and 
screened potentials, the Fourier coefficients v^ k „' k < are 
given by 
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respectively, where /5^ kn / k , are the Fourier coefficients of 
the pseudocharge density [see Eqs. (|A9[) - (j A13I) of Ap- 
pendix IA 2| . Note that for the unscreened potential, the 
term corresponding to q = (i.e., k = k' and G = 0) 
leads to a singularity which has to be considered carefully 
(details are given at the end of this section). 

The radial function u^kn'k' ^ s gi ven by [r< = 
min (r, r'), r> = max (r, r'), and R a is the radius of the 
atomic sphere] 

Ra 

vA(r) = J pA,(r')Gf(r,r')r' 2 dr' 
o 

+vA(R a )Pt(r), (20) 

where Gf is Eq. (|A7I) and Pi — r l / R l a for the unscreened 
potential or Gf is Eq. (fM|) and Pi = it (Xr) /it (XR a ) 
for the screened potential. In Eq. (|2U1) . 

VA (Ra) = 4^ E <ln^ Ta YL (5) 31 (|q| , 
G 

( 21 ) 

which is obtained by using the Rayleigh formula^ 

oo I 

j*' = 4tt J2 E ^ Cl<li r) YL (q) ^ (f ) (22) 

1=0 m=-t 

in the Fourier expansion of v^ k „/ k / [Eq. (HU], where jt 



is a spherical Bessel function^ 

E^ v is decomposed into its interstitial and atomic 
sphere parts: 
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where Q e = 1/R £ a and g"™ k "' k ' is Eq. (|AT5|) for the 
unscreened potential or Qg — (1/ie (XR a )) X / (2£ + 1)!! 
and g"„" k "' k ' is Eq. (|A17I) for the screened potential. 

As already mentioned above, the singularity which 
arises when q = [see Eq. (fT8|) ] needs to be consid- 
ered properly. Several methods to deal with this in- 
tegrable singularity when integrating into the Brillouin 
zone are available in the literatur o 27 ' 60 i 65 ~ — and have 
been used in very recent studiesi 61 ' 70 " — We adopted the 
simple scheme proposed by Spencer and Alavi£2> which 
consists of multiplying Eq. (H"8f by 1 — cos (| q| R c ), where 
R c = (3/ (4tt) N k tt) 1/3 with AT k being the number of k- 



points in the full Brillouin zone. In the real space this 
corresponds to multiplying Eq. (|5]) by the step function 
8(R C — |r — r'|)j££ By doing this, the term q = tends 
to a finite value: 



lim 

|qK0 



4tt 



(l-cos(|q| R c ))=2irR 2 c 



(27) 



which leads to a much more faster convergence (with re- 
spect to iV k ) of the integrations into the Brillouin zone. 

The screened potential has no singularity at q = [see 
Eq. (|19[) ]. nevertheless it is still useful to apply the same 
technique in order to accelerate further the convergence 
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of the integrations into the Brillouin zone. Multiplying 
Eq. ^ by the step function 9(R C — |r — r'|) means that 
in the reciprocal space Eq. (fP9|) should be multiplied by 

1 - e- XR < ^ A S in (|q| R c ) + cos (|q| R c ^j , (28) 
which becomes 1 — e~ XRc (XR C + 1) at q = 0. 



2. Valence-core and core-core terms 



By supposing that the core shells are closed (see Refs. 
l2~il and [23), the Legendre polynomial addition theorem^ 3 - 
can be used to simplify the calculation of the valence-core 
and core-core terms of the HF exchange energy [Eqs. ([3]) 
and The final expressions are given by 
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where H t {r,r') = (4tt/ (2£ + 1)) r^/ri^ 1 for the un- 
screened potential or ff«(r, r') = 47rAif (Ar<) fc^ (Ar>) for 
the screened potential. Cf 3 ™ 3 » „, are Gaunt coefficients 

[Eq. (TTT])] and £>"™ k/ were defined in Sec. H3 Note 
that in Eqs. (f2U| and (f3"U|). all integrations are inside the 
atomic spheres only, thus the cost for the calculation of 
these two terms is negligible compared to the valence- 
valence term. 



V™ = Wxct.w + «xi7,vc- F° r tri e present work we chose to 
implement the HF (and hybrid, see Sec. IIIip operator 
using a second variational procedure, which consists of 
using the semilocal (SL), LDA or GGA, orbitals as basis 
functions for the calculation of the matrix elements of 
the perturbation operator (ip^ h \v^ - Tne 
HF part is given by 



B. Hamiltonian 



The HF exchange operator for the valence orbitals is 
the sum of the valence-valence and valence-core terms: 
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(32) 



which are calculated using the same procedure as for the Eq. (1311) . The second variational procedure, which was 
HF exchange energy, but with /< kn „ k „ = Ck L *C" k » for 
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also adopted for the HF implementations in other LAPW 
code o 21 i 22 ' 27 i 61 leads to cheaper calculations, since in 
practice the number of orbitals ip^i^ which are used for 
the construction of the HF Hamiltonian matrix is much 
smaller than the number of LAPW basis functions. In 
the present implementation, the core electrons experience 
the semilocal potential, similarly as what is done in the 
fleur code, where the core electrons are taken from a 
previous semilocal calculation and kept frozen during the 
calculation with the hybrid functional^ 

III. SCREENED HYBRID FUNCTIONALS 

In screened hybrid functionals, the SR part of a frac- 
tion a x of semilocal exchange is replaced by SR HF 
exchanged 

£ xc = Eg + a x (£f - HF - £ X SR - SL ) , (33) 

where £ , ^ R " HF and E x R_SL are obtained by replacing the 
full (i.e., unscreened) Coulomb operator by the screened 
(i.e., SR) operator into the corresponding expressions. 
For the HSE functional^ the Coulomb operator was 
split into SR and LR components by using the error func- 
tion, however, for the present work we chose to split the 
Coulomb operator by using the exponential function: 

1 _ e - A l r - r 'l l-e- A l r - r 'l 
|r-r'| ~ |r-r'| + ^ |r-r'| ~ ' ' ' ' 

SR LR 

Figure Q] shows the SR and LR parts [Eq. ([3"I|] of the 
Coulomb potential 1/x = 1/ jr — r'|, and for compari- 
son, the same is shown when the error function is used 
to split 1/x [eric(px) = l-erl(px) is the complementary 
error function]. In both cases, the screening parame- 
ter is set to A = \i = 1. At x = 0, the values of the 
LR parts (l — e~ Xx ) jx and (1 — erfc(/xx)) fx are A and 
2[i 1 0F, respectively, thus these two ways of splitting the 
Coulomb operator lead to LR components which are not 
zero at x — 0. Sharper splitting schemes which lead to 
a LR component which is zero at x = consist of using, 
e.g., the erfgau function^ 6 - or simply the step function. 69 
We mention that for technical convenience, Shimazaki 
and Asai replaced e~ Xx by erfc((2/3)Aa:) in their pro- 
posed screened HF potentiali^Zr— Indeed, from Fig. [2] 
we can see that if A = (3/2) fi, the two splitting pro- 
cedures lead to very similar SR and LR parts. In this 
example, /i = 0.11 bohr -1 , which is the value used for 
the HSE06 functional^ 

In Eq. (J33|), E^ R - HF is given by Eqs. ©-CI]) with the 
Yukawa potential [Eq. ©] for v and E^ R ~ SL is given by 

£x SR - SL = -§(£) 1/3 £/> (r) 

xF x (s a (r))J(a a {r))d 3 r, (35) 




FIG. 1: (Color online) Plots of the SR and LR parts of the 
Coulomb operator 1/x, when split using the exponential (in 
blue) or error (in red) functions. 

where F x (s a ) [where s a — |Vp CT | / (2p tT kp ) with fcp = 

1 /3 

(67r 2 p CT ) ] is the enhancement factor of the semilo- 
cal exchange functional and J(a a ) [where a a — 
Aa/F x (s ct )/ (2fcp)] is a function, whose analytical form 
depends on the way the Coulomb operator is screened 
[J(a a ) = 1 for the unscreened operator]. In our case [SR 
part of Eq. ([3i])]. J(a a ) is given by§2, 

t, s , 2 2 8 1 

J{a a ) = 1 - -a a - -a a arctan — 
3 3 fl ff 

+ ^a 2 a (al + 3)ln^l + ^ . (36) 

Equation (|35|) is an approximation which was originally 
proposed by Iikura et oi.jSi but with the function J(a a ) 
for the error function. Recently, Akinaga and Ten-nc 82 
used Eq. (|35j) in conjunction with the Yukawa poten- 
tial as in the present work. (However, we note that in 
Refs. [8l] and[H, this is the LR part of the semilocal ex- 
change which was replaced by LR HF.) A more elegent 
way of calculating E x R ~ sh would be to use its expression 
in terms of the exchange hole (as done for HSE 39,83 ), 
and, for practical convenience, to find a mathematical 
form for the exchange hole such that an analytical inte- 
gration with the Yukawa potential is possible, as done in 
Ref. HH for the error function. This method has been 
used in Ref. [85| for the HSEsol functional, which is based 
on the PBEsol GGA functional^ We did not consider 
this possibility for the present work. 

For the semilocal terms in Eq. (|33p we have chosen 
PBE, 4 which is of the GGA form. In the following, this 
functional will be called YS-PBEO (where YS stands for 
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FIG. 2: (Color online) Plots of the SR and LR parts of the 
Coulomb operator 1/x, when split using the error (with fi = 
0.11 bohr -1 ) or the exponential [with A = (3/2) ^ = 0.165 
bohr -1 ] functions. 

Yukawa screened). In the literature, the unscreened ver- 
sion of this hybrid functional (recovered when A —¥ 0) is 
called PBE0) 36 i 37 for which the fraction of HF exchange is 
a x = 0.25 (see Ref. For A -> oo, YS-PBE0 reduces 
to PBE. The calculation of the total energy for hybrid 
functionals is given in Appendix [B] As already men- 
tioned in Sec. IIIB| the second variational procedure has 
been implemented, and the matrix elements of the per- 
turbation operator corresponding to Eq. ([55)) are given 
by 

(r n tK (e- HF - e- SL ) \rJi), (37) 

where the expression for wf^" SL = <5i? x R ~ SL /Sp a is given 
in Appendix [Cl 

IV. NUMERICAL RESULTS 

The calculations presented in this section were done 
with values for the parameters such that the results are 
well converged. The most important parameters are the 
number of k points for the integrations into the Brillouin 
zone, the size of the basis sets (first and second varia- 
tional procedures), G max and £ max in Eqs. (fT3l) and (1T41) . 
and < max in Eq. (fl"5|) . We will not discuss in detail the 
convergence of the results with respect to these parame- 
ters, but just mention the following: for the transition en- 
ergies and lattice constants, the number of orbitals used 
as basis functions for the second variational procedure is 
between two and six times larger than the number of va- 
lence bands in the system. The values of G max lie in the 



TABLE I: Total and exchange energies (in Ha) of He atom. 



wien2k Reference 



Functional 


— -Btot 




— -Etot 




LDAx" 


2.724 


0.853 


2.724 


0.853 


B88° 


2.863 


1.016 


2.863 


1.016 


PW91x a 


2.855 


1.005 


2.855 


1.005 


HF a 


2.862 


1.024 


2.862 


1.026 


HF 6 




0.998 




0.998 


HF C 




1.017 







"Obtained from exchange-only self-consistent calculations. The 
reference results are from Refs. [5^ and l9ll . 

''Evaluated with LDA (exchange and correlation) orbitals. The 
reference result is from Ref. 11 

Evaluated with B88PW91 orbitals. 



range 4—10 bohr -1 and for most calculations the value 
f max = 4 was used, which is more than enough most of 
the time. The size of the k-meshes will be mentioned 
below. 

We mention again that the computation of the HF 
Hamiltonian is very expensive, and for the systems we 
have considered this leads to computational times which 
are by one or two orders of magnitude larger than for 
semilocal functionals. Actually, the values of all param- 
eters mentioned above have a large impact on the com- 
putational time. 



A. Comparison with other codes 

1. HF energy 

As a first test of the correctness and accuracy of the im- 
plementation, we considered systems which do not con- 
tain core electrons, such that all electrons are treated 
self-consistently with the HF method. The He atom and 
solid LiH are two such systems for which highly accurate 
HF results are available in the literature. The LDAx 
(exchange-only LDA) orbitals were used as basis func- 
tions for the Hamiltonian of the second variational pro- 
cedure. 

The results for the He atom are shown in Table [I] The 
calculations were done in a fee cell with a lattice con- 
stant of 9.5 A which is large enough to make the inter- 
actions between the He atoms negligible. First, in order 
to have an idea of the accuracy that can be expected 
with wien2k, we did calculations with semilocal func- 
tionals (exchange only: LDAx, B88r& and PW913* 3 -) and 
compared them to accurate atomic results ! 89 ' 90 From the 
results we can see that an agreement at the mHa level 
can be reached, which is the target for the HF calcula- 
tions. The self-consistent HF results shown in Table U 
were obtained using 410 bands for the second variational 
procedure, which was enough to reach convergence and 
thus agreement with accurate atomic results^ However, 
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number of bands 

FIG. 3: (a) Total and (b) exchange energy of the He atom 
with respect to the values calculated with 410 bands. 



for the exchange energy E x , the agreement with the ref- 
erence result is not perfect. Actually, we can see in Fig. 
[3] that for a given number of bands, the error with re- 
spect to the (approximately) converged value is ten times 
larger for the exchange energy than for the total energy. 
For E toi , about 120 bands are necessary to reach con- 
vergence at the mHa level, while 410 bands are still not 
enough for E x (about 3000 LAPW basis functions are 
used for the first variational procedure) . It is known that 
the total energy converges faster than its components. In 
order to evaluate the effects due to self-consistency, the 
HF exchange energy was also evaluated using the LDA 
(with PW92 for correlation^) and B88PW91 (B88^ for 
exchange and PW91— for correlation) orbitals. From Ta- 
ble HI we can see that using LDA orbitals leads to an HF 
exchange energy whose magnitude is 26 mHa smaller, 
while using B88PW91 orbitals leads to a value which is 
much closer to the self-consistent one, which is maybe 
not surprising since the empirical parameter in B88 was 
determined by a fit to HF exchange energy of rare-gas 
atoms4^ This indicates that using the B88PW91 orbitals 
as basis functions for the second- variational Hamiltonian 
would be more efficient. 

In Ref. |74| (as well as in the Comment), well converged 
calculations on solid LiH (rocksalt structure) using Gaus- 
sian basis sets yield a value of E to t — —8.0645 Ha at the 
experimental geometry (4.084 A). Using a 6 x 6 x 6 k- 
mesh and 65 bands for the second variational procedure 
we obtained E tot — —8.0642 Ha. Increasing further the 
number of bands would lower the total energy and reduce 
the difference between the Gaussian and LAPW results. 
Therefore, as for the He atom, the HF energy calculated 
with the LAPW code agrees very well with the literature 



results. 



2. PBE0 calculations 

In Refs. [44] and HH calculations with the unscreened 
hybrid functional PBE0 were done within the projector 
augmented-wave ("VASP code) and LAPW (fleur code) 
methods, respectively. The implementation of the HF 
equations within the LAPW basis set as reported in Ref. 
[6X1 was done using another technique (mixed product ba- 
sis) as the one used in the present work (pseudocharge 
method). The integrations into the Brillouin zone were 
done with a 7 x 7 x 7 k-mesh for the semiconductors 
and insulators, while for the metals Li, Cu, and Rh a 
12 x 12 x 12 k-mesh was used. The results from Refs. |H 
and l6ll were done with a 12 x 12 x 12 k-mesh, however 
test calculations indicate that our results are converged 
within ~ 0.02 eV for the transition energies and ~ 0.002 
A for the lattice constants. 

Transition energies were calculated for six solids at the 
experimental lattice constant: Ar (fee, 5.260 A), C (di- 
amond, 3.567 A), Si (diamond, 5.430 A), GaAs (zinc 
blende, 5.648 A), MgO (rocksalt, 4.207 A), and NaCl 
(rocksalt, 5.595 A). The PBE0 results, as well as the PBE 
and experimental results, are given in Table ILT1 where we 
can see that the wien2k results agree very well with the 
fleur and VASP results. There are a few cases where 
the discrepancy is larger than 0.1 eV. For the r — » L 
transition in C, there is a difference of 0.12 eV between 
the wien2k and fleur values and in the case of NaCl, 
a disagreement of 0.15—0.2 eV with fleur is found for 
the r — > r and T — >• X transitions. Nevertheless, overall 
the agreement with the fleur and VASP codes for the 
PBE0 hybrid functional is clearly satisfactory, in partic- 
ular with VASP. Compared to the experimental values, 
the PBE0 functional clearly improves upon PBE, how- 
ever, some sizeable disagreements with experiment are 
still present, as for example for Ar and NaCl for which 
PBE0 underestimates the r — > T transition by about 
3 and 1.2 eV, respectively. In general, the tendency of 
the PBE0 functional is to overestimate small band gaps 
(e.g., GaAs) and to underestimate large band gaps (e.g., 
rare-gas solids) 4£ 

The lattice constant and bulk modulus of a few se- 
lected compounds, namely, Li (bec), C (diamond), Si 
(diamond), Cu (fee), Rh (fee), LiF (rocksalt), BN (zinc 
blende), and SiC (zinc blende) were calculated using the 
PBE and PBE0 functionals. The results are shown in Ta- 
ble Hill together with the values obtained with the VASP 
code^ and the experimental data, which were corrected 
for the zero-point anharmonic expansion. 85 By compar- 
ing the wien2k and VASP results, we can see that ex- 
cellent agreement between the two codes are obtained 
both for the PBE and PBE0 functionals. The largest 
discrepancy in the lattice constant is found for Si, where 
a difference of 0.007-0.01 A is found for PBE and PBE0. 
From Table Hill we can see that there is also a good agree- 
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TABLE II: Transition energies (in eV) obtained with the PBE, PBEO, and YS-PBEO (A = 0.165 bohr -1 ) functional. 



Solid 


Transition 




WIEN2K 






VASP" 




FLEUR 6 


Expt. c 


PBE 


PBEO 


YS-PBEO 


PBE 


PBEO 


HSE06 


PBE 


PBEO 


Ar 


r - 


+ r 


8.69 


11.09 


10.36 


8.68 


11.09 


10.34 


8.71 


11.15 


14.2 


C 


r - 


->■ r 


5.59 


7.69 


6.94 


5.59 


7.69 


6.97 


5.64 


7.74 


7.3 




r - 


+ X 


4.76 


6.64 


5.91 


4.76 


6.66 


5.91 


4.79 


6.69 






r - 


+ L 


8.46 


10.76 


9.97 


8.46 


10.77 


10.02 


8.58 


10.88 




Si 


r - 


-> r 


2.56 


3.95 


3.30 


2.57 


3.97 


3.32 


2.56 


3.96 


3.4 




r - 


+ X 


0.71 


1.91 


1.31 


0.71 


1.93 


1.29 


0.71 


1.93 






r - 


4- L 


1.53 


2.86 


2.23 


1.54 


2.88 


2.24 


1.54 


2.87 


2.4 


GaAs 


r - 


-> r 


0.53 


1.99 


1.39 


0.56 


2.01 


1.45 


0.55 


2.02 


1.63 




r - 


+ X 


1.46 


2.66 


2.08 


1.46 


2.67 


2.02 


1.47 


2.69 


2.18, 2.01 




r - 


+ L 


1.01 


2.35 


1.74 


1.02 


2.37 


1.76 


1.02 


2.38 


1.84, 1.85 


MgO 


r- 


+ r 


4.79 


7.23 


6.49 


4.75 


7.24 


6.50 


4.84 


7.31 


7.7 




r - 


+ X 


9.16 


11.58 


10.83 


9.15 


11.67 


10.92 


9.15 


11.63 






r - 


+ L 


7.95 


10.43 


9.68 


7.91 


10.38 


9.64 


8.01 


10.51 




NaCl 


r - 


+ r 


5.22 


7.29 


6.61 


5.20 


7.26 


6.55 


5.08 


7.13 


8.5 




r - 


->■ X 


7.59 


9.80 


9.06 


7.60 


9.66 


8.95 


7.39 


9.59 






r - 


L 


7.33 


9.40 


8.70 


7.32 


9.41 


8.67 


7.29 


9.33 





"Reference 44 (see Erratum for HSE06 results). 
'Reference 6l]. 

' The references for the experimental values are given in Table 1 
of Ref. IU. 



TABLE III: Equilibrium lattice constants a (in A) and bulk moduli B (in GPa) obtained with the PBE, PBEO, and YS-PBEO 
(A = 0.165 bohr -1 ) functionals. The experimental values, which are corrected for the zero-point anharmonic expansion, are 
from Ref. jj. 

WIEN2K vasp° 



PBE PBEO YS-PBEO PBE PBEO HSE06 Expt. 



Solid 


do 


Bo 


Oo 


Bo 


ao 


Bo 


ao 


Bo 


ao 


Bo 


ao 


Bo 


ao 


Bo 


Li 


3.434 


13.9 


3.464 


13.1 


3.467 


12.6 


3.438 


13.7 


3.463 


13.7 


3.460 


13.6 


3.453 


13.9 


C 


3.575 


435 


3.549 


475 


3.554 


467 


3.574 


431 


3.549 


467 


3.549 


467 


3.553 


455 


Si 


5.476 


89.0 


5.443 


99.4 


5.459 


96.5 


5.469 


87.8 


5.433 


99.0 


5.435 


97.7 


5.421 


101 


Cu 


3.631 


141 


3.630 


131 


3.654 


119 


3.635 


136 


3.636 


130 


3.638 


126 


3.595 


145 


Rh 


3.830 


256 


3.787 


292 


3.799 


280 


3.830 


254 


3.785 


291 


3.783 


288 


3.794 


272 


LiF 


4.069 


67.2 


4.008 


70.2 


4.035 


67.3 


4.068 


67.3 


4.011 


72.8 


4.018 


72.7 


3.972 


76.3 


BN 


3.628 


374 


3.601 


407 


3.607 


401 


3.626 


370 


3.600 


402 


3.600 


402 


3.592 


410 


SiC 


4.384 


213 


4.352 


242 


4.361 


236 


4.380 


210 


4.347 


231 


4.348 


230 


4.346 


229 



"Reference (see Erratum for HSE06 results). 



ment between the two codes for the bulk modulus. On 
average, the hybrid functional PBEO improves over the 
GGA PBE for the lattice constant and bulk modulus of 
semiconductors and metals as shown in Ref. [lil 



3. YS-PBEO calculations 

As mentioned in Sec. Mil (see Fig. [5]), choosing 
A = (3/2) /i in Eq. (13"4"|) leads to a splitting of the 
Coulomb operator which is very similar to the one ob- 



tained by using the error function with a given /ziiZr— 
In the HSE06 functional^ /x is fixed to 0.11 bohr -1 and 
in order to see how well the YS-PBEO functional can 
reproduce the HSE06 transition energies (see Erratum 
of Ref. S3), calculations with A = (3/2)0.11 = 0.165 
bohr -1 were done. From the results shown in Table ILTl 
we can see that the agreement between HSE06 (vasp) 
and YS-PBEO is as good as it was for PBEO with dif- 
ferences smaller than 0.03 eV in most cases. Compared 
to PBEO, the screened hybrid functionals lead to bet- 
ter (worse) agreement with experiment for small (large) 



TABLE IV: Band gap (in eV) and Cu EFG (in 10 21 V/m 2 ) 
of CU2O calculated at the experimental lattice constant 
(4.27 A). 

EFG 



ivietnoa 




Band gap 


TV-it q1 
±Ot Ell 


p-p 


d-d 




- 


IjJJA 




U.OO 


— O.O 


— 10.U 


1U.0 










O K Q 


—0.0 


1 a a 
— 10.4 


1U.0 




- 
- 


DOorWyl 




U.00 


— 0.O 


1 d A 

— 10.4 


1U.0 






hvyor VV91 




0.0 < 


—0.0 


1 T A 

— 1 ( .4 


10.D 




- 
- 


T T\ A 1 T T /T?T T T T 

LDA+f (r_L_L, U — 


4 ev J 


0.00 


—0.1 


1 a 1 
— 10. 1 


9.0 






rr>A 1 rr /rriT t rr 


ev ) 


U.oU 


— 0.0 


1 a a 
— 10.4 


y.o 






T T\ A 1 rr /PT T TT 

LDA+tv (r-L-L, L7 = 


12 ev) 


0.91 


—7.6 


— 16.5 


8.8 


O 


■: 


LiJA+t/ ^Alvlr , - 


— 4 ev J 


U.Oo 


— 4.0 


— 10. U 


11 n 

IX. u 






T T\ A 1 f 7" / A TV .f T7 1 I" T 

LJJA+(V (A Mr , U = 


= e V J 


0. (9 


— ^.0 


ICQ 
— lO.Z 


1 0.4 




■: 


LDA+f/ (AMF, C7 = 


= 12 eV) 


0.94 


0.6 


-16.6 


17.0 




— 


PBEO (onsite) 




0.79 


-3.4 


-16.4 


12.8 






PBEO 




2.77 


-8.5 


-19.5 


10.8 






YS-PBEO 




1.99 


-8.3 


-19.3 


10.8 






pseudo-SIC a 




1.80 












B3LYP b 




2.1 












HSE (a x = 0.275) c 




2.12 












scGW d 




1.97 












Expt. 




2.17 e 


9.8 / 











a Reference 9i|. 
'Reference 9a. 
c Reference Hoi 
d References 
e References 

^Only the magnitude 
Q 106.107 



102 



105 



is known. Calculated using Q ( 63 Cu) 




FIG. 4: (Color online) Density of states of CU2O calculated 
with different functionals. The Fermi energy is set at E — 
eV. 



band gaps (see also Ref. |47| ). 

The YS-PBEO results for the lattice constant and bulk 
modulus are shown in Table Hill The agreement between 
the HSE06 and YS-PBEO results is fairly good in cases 
like Li or C, while larger differences can be seen for Si 
(0.024 A), LiF (0.017 A), Cu (0.016 A), and Rh (0.016 
A) . An important contribution to these differences in the 
lattice constant between the HSE06 and YS-PBEO values 
could be attributed to the different schemes used for the 
screening of the semilocal exchange term [Eq. (|35|) ]. For 
YS-PBEO, the method of Iikura et alM is used, while 
in HSE06 the screened exchange energy is obtained by 
integrating a model of the exchange hole . 39 ' 83 However, 
it seems that using one of the scheme or the other has 
very little influence on the transition energies as shown 
above. 



B. Cu 2 

Cuprous oxide (CU2O) is a semiconductor which has 
been used in many applications (e.g., catalysis and photo- 
voltaics). Its structure is cubic (space group Pn3m) and 
the unit cell, which has a lattice constant of 4.27 A, 93 



contains six atoms. In this structure, shown in Fig. 1 
of Ref. [94|, the O atoms are fourfold coordinated by Cu 
atoms, whereas the Cu atoms are linearly coordinated by 
O atoms. Formally Cu has a valency of +1 and the Cu- 
3d shell in CU2O is full, therefore the correlation effects 
in the Cu-3<i shell should not play an important role as 
it is the case for CuO^S- 

Many experimental and theoretical studies on CU2O 
have been done. On the theoretical side it has been 
shown that the semilocal approximations underestimate 
the band gap as expected (see Refs. [9fj| and [93 for col- 
lections of previously done calculations), but also the Cu 
electric-field gradient (EFG), 98 which is a ground-state 
property derived from the electron density. LDA+J7 
(or GGA+U) improves only slightly over the semilo- 
cal approximations , 97 ' 98 while the pseudo self-interaction 
method (pseudo-SIC), 99 the hybrid functionalsffi 100 ' 101 
and self-consistent GW (scGW)^- provide band gaps in 
much better agreement with experiment. Actually, the 
results for the EFG show that the semilocal and LDA+C7 
methods do not provide an accurate description of the 
occupied states. 

In Table IIV1 we show the results for the band gap 
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and EFG obtained with the hybrid junctionals PBEO 
and YS-PBEO (A = 0.165 bohr" 1 ), which were obtained 
with a mesh of 5 x 5 x 5 k-points. The calculations 
with the semilocal (LDA^ PBE,i BSSPWOl^ and 
EV93PW912JH), LDA+U [fully localized limit (FLL)±2i 
and around mean- field (AMF) 104 versions], and onsite 
PBE0 5; methods (results in Table [TV) were done with 
a 12 x 12 x 12 k-mesh. The radii of the Cu and O 
atomic spheres are 1.84 and 1.63 bohr, respectively. 
LDA, PBE, and B88PW91 give values for the band gap 
(~ 0.5 eV) and EFG (~ -5.5 x 10 21 V/m 2 ) which 
are much smaller than the experimental values (above 
2 eV for the band gap^ and 9.8 x 10 21 V/m 2 for the 
EFG 106 i 107 ). EV93PW91 improves for the EFG with a 
value of —6.6 x 10 21 V/m 2 , but not for the band gap con- 
trary to what was reported for many other solids in Ref. 
1 1081 . LDA+J7 slightly improves the results for the band 
gap and its two versions, FLL and AMF, lead to the same 
value for a given value of the Coulomb parameter U (the 
exchange parameter J has been fixed to 0.95 eV). How- 
ever, this improvement is minor and even with U = 12 eV 
the band gap remains well below the experimental value. 
For the EFG, FLL and AMF lead to different trends. An 
increase of U leads to an increase of the magnitude of 
the EFG with FLL, while the opposite is obtained with 
AMF, which yields a positive value for U = 12 eV. In 
Table IIV1 the p-p and d-d components (inside the Cu 
atomic sphere) of the EFG are also shown. As expected, 
the change in the EFG due to U comes mainly from the 
d-d part. The onsite PBEO method slightly improves the 
results for the band gap (0.8 eV), but significantly de- 
creases the EFG (-3.4 x 10 21 V/m 2 ). Overall, the FLL 
version of LDA+U leads to a moderate improvement over 
the semilocal functionals, while AMF and onsite PBEO 
behave similarly by reducing the magnitude of the EFG. 

The results obtained with PBEO and YS-PBEO are in 
much better agreement with experiment. In particular, 
the screened YS-PBEO functional leads to a band gap 
of 1.99 eV, which is very close to the experimental value, 
and an EFG of —8.3 x 10 21 V/m 2 which is much closer to 
experiment compared to the values obtained with other 
functionals. PBEO leads to a band gap which seems to 
be too high and an EFG very similar to YS-PBEO. Fo- 
cusing now on the PBE and PBEO results for the EFG, 
we can see from the decomposition of the EFG (Table 
IIV[) that the increase in magnitude of the EFG by going 
from PBE to PBEO comes mainly from the p-p compo- 
nent. By decomposing further the p-p component, we 
could see that the sub-component from the Cu-4p states 
(which actually originate mainly from a re-expansion of 
the 0-2p tails) is more negative than the total p-p and 
that the sub-component from the low-lying (~ —5 Ry) 
semicore Cu-3p states is small and positive. From this 
we could also determine that the more negative PBEO 
p-p component come half and half from the Cu-3p and 
Cu-4p states. 

Figure @] shows the density of states (DOS) of CU2O 
for a few selected functionals. We can see that in the 



energy range between —8 and — 5 eV below the Fermi 
energy (set at E = eV), most of the DOS comes from 
0-2p electrons. The DOS between —4 and eV is en- 
tirely due to Cu-3g? states, while above the band gap, the 
different partial DOSs actually represent Cu-4s states. 
By comparing the different functionals, we can observe 
that the 0-2p peaks are higher in energy (closer to the 
Cu-3d states) for the FLL version of LDA+C7. Also, the 
LDA+U and hybrid methods shift the main Cu-3d peaks 
down in energy. 

Compared to the other hybrid results from the litera- 
ture (also shown in Table HV)) . we can see that the HSE 
(with a x = 0.275) band gap of 2.12 eV-lSi is close to the 
YS-PBEO value of 1.99 eV, as expected from the results 
obtained in Sec. II V A 31 The B3LYP band gap of 2.1 eV 
reported in Ref. |96j is much smaller than our PBEO value 
of 2.77 eV. This is mainly due to the smaller amount of 
HF exchange a x in B3LYP (0.2 for B3LYP versus 0.25 
for PBEO). 



V. SUMMARY 

The implementation of unscreened and screened hy- 
brid functionals into the wien2k code, which is based on 
the LAPW basis set, has been presented. The screening 
is based on the Yukawa potential for which the expansion 
in spherical harmonics has a simple expression. Also, it 
was possible to calculate analytically all integrals which 
were necessary for the derivation of the various formulas 
for the pseudocharge method. In order to check the va- 
lidity of the implementation, first, test calculations were 
done on systems which do not contain core electrons, 
such that the total Hartree-Fock energy could be com- 
pared with benchmark results from the literature. As a 
further test of the implementation, the band gap and lat- 
tice constant of several solids have been calculated with 
the hybrid functionals PBEO and YS-PBEO which are 
based on the GGA functional PBE. The results are in 
very good agreement with the results obtained by other 
codes. Noticeably, for the screened hybrid functional YS- 
PBEO, it was possible to find a value of the screening 
parameter A such that the results are very close to the 
results of HSE06, whose screening is based on the error 
function. Finally, we applied the hybrid functionals to 
the semiconductor CU2O. The results obtained with the 
unscreened PBEO and screened YS-PBEO for the band 
gap and EFG are much more accurate than the results 
obtained with semilocal and LDA+U functionals. 



Appendix A: Pseudocharge method 

In this appendix, the basic formulas of the pseu- 
docharge method are given (Sec. IA 1|) . as well as explicit 
expressions used for the HF energy (Sec. I A 21) . 
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Basic formulas 



or 



In all-electron calculations, the charge density p has 
large oscillations near the nuclei, therefore its Fourier ex- 
pansion will converge slowly, making the calculation of 
the potential v generated by p with Fourier transforms 
inefficient. The idea of the pseudocharge method^ 8 - is to 
replace the charge density p inside the atomic spheres S a 
by a smoother one (p) such that the Fourier expansion of 
P = Ppw + P converges faster, ppw is the continuation in- 
side the spheres of the plane waves (PW) representation 
of the charge density and p = ^ Q ~p a is zero in the inter- 
stitial region. Such a scheme is possible since the poten- 
tial in the interstitial region created by the charge inside 
the spheres depends only on the multipole moments qf m 
which are defined as follows for the unscreened [Eq. ©] 
and screened [Eq. ©] potentials: 



Y; m (r)r e p(r)d 3 r, 



Gf(r,r') = 47rAi*(Ar<)fc(Ar>) 

for the unscreened and screened potentials, respectively. 
dG a /dn' is the normal derivative of G a at the sphere 
boundary. 



2. Explicit expressions for the Hartree-Fock energy 



In Eqs. (|18[) and (fl9)). the Fourier coefficients of the 
pseudocharge density are given by 



X<"1 _ crG , -<rq 

f'nkn'k' — f'nkn'k' T P„kn'k' 



(A9) 



(Al) 

where (unscreened case) 



Qlm = (2 ^/ )!! / YL (f ) H (Ar) P (r) d 3 r. (A2) 

s Q 

Therefore, ~p a should be chosen such that inside the 
spheres, the multipole moments of p are equal to the 
multipole moments of the true charge density p. After 
having determined p, the potential in the interstital re- 
gion is calculated with 



G M 



(A3) 



and 



^nkn'k' ^ 



cell 

EE 



4tt^;^ (2^ + 2p + 3)!! (-i) e 



Rl + p+l {2e+1) 



, jl+p+1 (|q| Raj „~iq-T aV /-\-aaukn'k' 

■ p+ i e i£ m (q)q em 



(A10) 



or (screened case) 



_ CTq = 47r yy A ' +p+1 H) 

r nkn'k' q / y / y 



n ^f^n +p+1 (XR a )(2£ + i)\\ 

, jt+p+1 (|q| Rot) -iq-T aV f-~r\—aankn'k' 

■ p+1 e ^m(,qjg£ m 



(r)=47r^ 



G 



— ^ e iG * 

G| 2 + A 2 



(A4) 



for the unscreened and screened potentials, respectively. 
Then, inside the atomic sphere S Q , the potential is the 
solution of a Green function problem: 

v a (r) = J p(r')G°(r,r')dV 

So 

* u, (r') ^ (r,r')sin#Wd0',(A5) 



where the Green function is given bj™ 

CO £ 



58,109 



(r, r') = E E G " ( r ' r ') ^/m(f ')^/m(f ), (A6) 

where 



t=0 m=-l 



4tt ri- f rl e+1 



2e + i r { +1 \ Rl 



2£+l 



(A7) 



(All) 



In Eqs. (|A10[) and (|A1 1|) . p is an integer which is chosen 
such that I + p is fixed. 58 For q = 0, Eqs. (|A10J) and 
(jAllj) reduce to 



cell 

PnWk = ^E«oo nkn ' k (A12) 



and 



r nkn'k o / j 



n ^ {2p + 3y.H p+1 {\R a ) 



— atrnkn'k 

</oo ) 



respectively. In Eqs. (|A10|) - (jA13|) . 



(A13) 



—acrnkn'k' acrnkn'k' PW, acrnkn'k' / * , \ 

<«m — Win — <«m > 

where q1^ nkn k and q^ ,aankn k are the multipole mo- 
ments of /Q^kn'k' inside the spheres and of the contin- 
uation of the PW representation of p„ k „/ k / inside the 
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spheres, respectively, whose expressions are given by xe lq ' T °y/ m (q) Pnkn'k' (A16) 

9"m" k "' k ' = E E T a/nkn>k' m for tlie unscrcene d potential and 

R oankn'k' _ \ Ll + \ \ rjihh&Mm 

/"lm \l / j / j ± aanVri'V! 

U h^ u h^ r ) ri+2(ir ( A15 ) w>/iJ> 

o Kq 



x / U £Vr)u%UrMAr)r 2 *(A17) 



9«« 



_ ^47riX +2 ^+i(|q|^) 



E 

G 



4 Wiarakn ' k ' = E (111 ^) ^-i ( Ai? «) - H h-i (|q| ^) (Aiia)] 4 ^ f a( f + 1 f e iq - T °y; m (q) p3U (ais) 



for the s creen ed pot entia l. When k = k', the term G = + a x (E^ F - E^ L ) , 

in Eqs. (|A16I) and (|A18|) reduces to 

feVS^rC'k (A19) 

3 where T s is the kinetic energy of the electrons and 



and 

?2 



^ ^[^ PnL^ (A20) ,c oul (r)= / ^1.^ - ^ _^_ | (B2) 



respectively. 

Appendix B: Total energy 



crystal 
— * 

crystal 



crystal 



|— — — — d r — 53 | Tq _ r | ' C^) 

tional, the total energy is given by (spin-unpolarized crystal 



form) 



■y r ^ ccl1 are ^ ne Coulomb and Madelung potentials, respectively, 

^tot =Ts+ 2 / l 'Coui(r)p(r)(i 3 r--^Z a u;(T a ) By using the sum of the eigenvalues 



cell 



n c ,£ c ,m. 



53 £ n c 4m c +E w »k£nk = T s + / u Cou i(r)p(r)d 3 r+ / u^(r)p(r)d 3 r+Q; x 2£££ v + £^ c - / ^ L (r)p va i(r)d 3 r 

™' k cell cell \ cell 

! 

where p va i is the valence electron density (the core elec- can be rewritten as 
trons experience the semilocal potential), the total energy 



E to t = 5Z e «cfcm c + 5Z w "k£nk ~ ^ v Cou\(r)p{r)d 3 r - -^2z a v^(T a ) ~ v^(r)p(r)d 3 



r 



"> k cell " cell 
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+E- 



SL 

xc 



R 



HF 



E. 



HF 



+ / v^{v)p val {r)d A r-E : 



SL 



(B5) 



cell 



For a screened hybrid functional, the exchange-only 
terms are simply replaced by their SR counterparts. The 
use of the second variational procedure allows us to write 
the sum of the valence eigenvalues in the following way: 



n . k 



n,k 



|2 SL 
^mk 



(2K 



HF 



V-2 = 



36 x 



2 A 



1 A 



1 A V 2 / 



36 S 



2 b 
1 A t 
2b* p 3 
3 A si 



2 6 3 p 5 / J 



-44 uy 2 



6 S 4 6 4 p 3 



,dH x 
ds 



El 



F 

1 X 



1/2 ' 



F 1/2 H 



(C3) 



HF 



^ L (r)p val (r)d 3 r , (B6) 



cell 



where c™ k are the coefficients of the expansion of ^„k 
(V>nk = EmCkV4k)- From E<1- dESJ), the valence- 
valence HF exchange energy E^ v can be calculated, thus 
avoiding the use of Eq. ([2]) , which is the most expensive 
component of the total energy to calculate. 



Appendix C: Functional derivative of 

The functional derivative of £f R ~ SL [Eq. (J35J)] for the 
spin-unpolarized case is given by 



t>3 



1A 2 s 2 
6 b 2 pVs 



1 A 2 



1 A 2 st 



3 6 2 flVa 4 „io/3 



where 6 = 2 (3tt 2 ) 1/3 , s = |Vp| / (2 ( 
Vp ■ V |Vp| , and iJ x = (1/s) cLF x /ds. 



(C4) 

3^ V V /3 U = 



,,SR-SL 



3 /3 



4 V 7T 



1/3 



UlJ + U2^- + U 3TT , (CI) 

da da z I 
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1_ 3 P x b 2 ^^^^^^ 5 6 3 p8/ 3 J ds 
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